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Abstract 

In  recent  work  [Phys.  Rev.  E  68  (2003)  025103],  it  was  shown  that  the  requirement  of  Galilean  invariance  determined 
the  form  of  the  H  function  used  in  entropic  lattice  Boltzmann  models  for  the  incompressible  Navier-Stokes  equations  in  D 
dimensions.  The  form  obtained  was  that  of  the  Burg  entropy  for  D  —  2,  and  the  Tsallis  entropy  with  q—  1  —  2/D  for  D^2. 
The  conclusions  obtained  in  that  work  were  restricted  to  particles  of  a  single-mass  and  speed  on  a  Bravais  lattice.  In  this  work, 
we  generalize  the  construction  of  such  Galilean-invariant  entropic  lattice  Boltzmann  models  by  allowing  for  certain  models 
with  multiple  masses  and  speeds.  We  show  that  the  required  H  function  for  these  models  must  be  determined  by  solving  a 
certain  functional  differential  equation.  Remarkably,  the  solutions  to  this  equation  also  have  the  form  of  the  Tsallis  entropy, 
where  q  is  determined  by  the  solution  to  a  certain  transcendental  equation,  involving  the  dimension  and  symmetry  properties 
of  the  lattice,  as  well  as  the  masses  and  speeds  of  the  particles. 

©  2004  Published  by  Elsevier  B.V. 
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1.  Introduction 


Lattice  Boltzmann  models  of  fluids  [2,3]  evolve  single-particle  distribution  functions  in  discrete-time  steps  on  a 
regular  spatial  lattice.  That  is,  velocity  space  is  discrete,  and  is  comprised  of  (possibly  linear  combinations  of)  the 
lattice  vectors  themselves.  In  spite  of  this  very  radical  simplification  of  the  Boltzmann  equation,  it  has  been  shown 
that  the  incompressible  Navier-Stokes  equations  emerge  unscathed  in  the  limit  of  small  Mach  and  Knudsen  numbers. 

In  the  most  general  situation,  the  collection  of  velocities  that  is  retained  do  not  all  have  the  same  magnitude,  and  we 
denote  them  by  caj,  where  the  index  a  is  associated  with  the  magnitude,  ca  =  |caj|,  and  j  enumerates  the  velocity 
vectors  with  that  speed;  velocities  with  the  same  index  a  are  said  to  be  in  the  same  speed  class.  The  single-particle 
distribution  corresponding  to  lattice  vector  c aj  at  lattice  position  x  and  time  step  t  is  denoted  by  Naj(\,  t).  The 
simplest  lattice  Boltzmann  models  employ  a  collision  operator  of  BGK  form  [4],  so  that  their  evolution  equation  is 


N, 


raj(x  +  C aj,  t  +  At)  =  NaJ(x,  t )  +  t )  -  NaJ(x,  /)] 


(1) 


for  j  =  1, . . .  ,  ba,  for  each  speed  a.  Here  ba  is  the  number  of  velocities  with  speed  ca,  N^ef  (x,  t )  is  a  specified 
equilibrium  distribution  function  that  depends  only  on  the  values  of  the  conserved  quantities  at  a  site,  and  r  is 
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a  characteristic  collisional  relaxation  time.  Using  a  discrete-velocity  version  of  the  Chapman-Enskog  analysis 
[3],  we  shall  show  that  the  mass  and  momentum  moments  of  the  distribution  function  obey  the  incompressible 
Navier-Stokes  equations  for  certain  choices  of  equilibrium  distribution. 

The  viscosity  that  appears  in  the  Navier-Stokes  equations  obtained  from  these  models  is  proportional  to  r  -  1  /2. 
To  lower  viscosity  and  thereby  increase  Reynolds  number,  practitioners  often  “over-relax”  the  collision  operator 
by  using  values  of  x  in  the  range  (1/2,  1].  Though  the  method  is  guaranteed  to  be  numerically  stable  for  r  >  1,  no 
such  guarantees  apply  when  x  <  1,  and  the  method  is  fraught  with  numerical  instabilities,  which  limit  the  highest 
Reynolds  numbers  attainable. 

In  an  effort  to  understand  and  thereby  avoid  these  instabilities,  there  has  been  much  recent  interest  in  entropic 
lattice  Boltzmann  models  [5-7].  These  models  are  motivated  by  the  fact  that  the  loss  of  numerical  stability  is  due  to 
the  absence  of  an  //-theorem  [7].  Numerical  instabilities  evolve  in  ways  that  would  be  precluded  by  the  existence  of 
a  well-behaved  Lyapunov  function.  The  idea  behind  entropic  lattice  Boltzmann  models  is  to  specify  an  H  function, 
rather  than  just  the  form  of  the  equilibrium  distribution;  of  course,  the  latter  is  that  which  extremizes  the  former. 
The  evolution  will  be  required  never  to  decrease  H,  yielding  a  discrete-time  //-theorem;  this  is  to  be  distinguished 
from  other  discrete  models  of  fluid  dynamics  for  which  an  //-theorem  may  be  demonstrated  only  in  the  limit  of 
vanishing  time  step  [8],  or  not  at  all. 

To  ensure  that  collisions  never  decrease  H,  the  characteristic  collision  time  r  is  made  a  function  of  the  incoming 
state  by  solving  for  the  smallest  value  xm\n  <  1  that  does  not  increase  H.  The  value  then  used  is  r  —  rmj n//c,  where 
0  <  k  <  1.  It  has  been  shown  that  the  expression  for  the  viscosity  obtained  by  the  Chapman-Enskog  analysis  will 
approach  zero  as  k  approaches  unity  [5-7],  Thus,  the  entropic  lattice  Boltzmann  methodology  allows  for  arbitrarily 
low  viscosity  together  with  a  rigorous  discrete-time  //-theorem,  and  thus  absolute  stability.  The  upper  limit  to  the 
Reynolds  numbers  attainable  by  the  model  is  therefore  determined  by  loss  of  resolution  of  the  smallest  eddies, 
rather  than  by  loss  of  stability  [7,9-1 1]. 

In  an  earlier  paper,  we  constructed  entropic  lattice  Boltzmann  models  for  the  incompressible  Navier-Stokes 
equations  that  are  Galilean-invariant  to  second  order  in  Mach  number,  and  we  showed  that  the  requirement  of 
Galilean  invariance  makes  the  choice  of  H  function  unique.  More  specifically,  we  showed  that  the  required  function 
has  the  form  of  the  Burg  entropy  [12]  in  two  dimensions,  and  the  Tsallis  entropy  in  higher  dimensions.  These 
conclusions  were  based  on  single-speed  lattice  Boltzmann  models  on  Bravais  lattices. 

In  this  work,  we  generalize  the  construction  of  such  Galilean-invariant  entropic  lattice  Boltzmann  models  by 
allowing  for  the  treatment  of  certain  models  with  multiple  particle  masses  and  speeds.  We  show  that  the  required 
H  function  for  these  more  general  models  must  be  determined  by  solving  a  certain  functional  differential  equation. 
Remarkably,  the  solutions  to  this  equation  also  have  the  form  of  the  Tsallis  entropy,  where  q  is  determined  by  the 
solution  to  a  certain  transcendental  equation.  This  equation  involves  the  dimension  and  symmetry  properties  of  the 
lattice,  as  well  as  the  masses  and  speeds  of  the  particles  in  the  model. 

In  Section  2,  we  describe  the  lattices  used  for  our  multispeed  lattice  Boltzmann  models.  In  particular,  we  introduce 
notation  for  sums  of  outer  products  of  velocity  vectors,  which  play  an  important  role  in  the  Mach  number  expansion 
of  the  equilibrium  distribution  and  the  Chapman-Enskog  analysis.  To  make  the  analysis  tractable,  we  restrict  our 
attention  to  multi-speed  models  for  which  each  speed  class  separately  satisfies  the  isotropy  requirements.  In  Section  3, 
we  specify  the  form  of  the  equilibrium  distribution  function,  assuming  only  that  the  H  function  is  of  trace  form, 
and  we  derive  the  Mach  number  expansion  of  this  equilibrium.  In  Appendix  A  we  construct  a  lattice  BGK  kinetic 
equation  with  this  equilibrium  distribution  and  apply  the  Chapman-Enskog  analysis  to  it,  solving  it  in  the  asymptotic 
limit  of  small  Knudsen  and  Mach  number.  The  resulting  hydrodynamic  equations  are  presented  in  Section  4.  These 
are  of  Navier-Stokes  form,  and  the  derivation  yields  the  equation  of  state  for  the  pressure,  the  viscosity,  and  the  g 
factor  that  may  multiply  the  convective  derivative.  In  Section  5,  we  examine  the  requirement  that  g  =  1,  needed  to 
restore  Galilean  invariance  in  the  context  of  four  examples.  The  first  is  the  single-speed  model,  in  order  to  show  that 
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our  methodology  is  able  to  recover  the  previously  known  results;  the  second  is  a  single-speed  model  augmented  by 
the  addition  of  rest  particles  of  the  same  mass;  the  third  allows  the  rest  particles  to  have  a  different  mass  front  the 
moving  particles;  and  the  fourth  is  the  general  case,  subject  to  the  above-described  restriction  on  our  lattice. 


2.  Description  of  lattice 


As  noted  above,  velocity  vectors  are  grouped  into  “speed  classes”  based  on  their  magnitudes.  We  may  associate 
these  magnitudes  with  speeds,  since  the  unit  of  time  is  taken  to  be  At  =  1.  We  denote  the  velocity  vectors  by 
c aj,  and  assume  that  these  are  (linear  combinations  of)  lattice  vectors.  Here  a  is  an  index  denoting  a  certain  speed 
class,  and  j  indexes  the  vectors  within  that  class.  All  lattice  vectors  within  a  speed  class  have  the  same  magnitude 
|ca,y|  =  ca,  and  are  associated  with  particles  of  the  same  mass,  ma. 

Sums  of  outer  products  of  the  velocity  vectors  arise  frequently  in  the  analysis  of  lattice  Boltzmann  models.  Within 
a  speed  class,  these  sums  are  denoted  by 

P; '  n 

b«,a  =  (2) 

j 

so  that,  in  particular,  bo,«  =  ba.  We  assume  that  these  quantities  vanish  for  odd  n,  and  that  they  are  isotropic  tensors 
for  even  n  <4.  That  is,  we  assume  that 

Ha  =  Hal,  (3) 


but  we  shall  not  specialize  to  this  case  in  what  follows.  The  assumption  that  the  outer  products  of  lattice  vectors  for 
even  n  <  4  be  isotropic  tensors,  separately  within  each  speed  class,  is  restrictive.  It  eliminates  from  our  consideration 
certain  interesting  models  (such  as  the  D2Q9,  D3Q15  and  D3Q1 9  models  [3])  in  which  these  tensors  are  not  isotropic 
within  each  speed  class,  but  are  so  when  combined  in  weighted  sums  across  speed  classes.  Nevertheless,  it  admits 
an  important  class  of  multi-speed  models  while  simplifying  the  analysis  considerably,  so  we  restrict  attention  to 
that  case  in  this  paper. 


3.  Equilibrium  distribution 

The  conserved  quantities  that  we  shall  consider  are  the  mass  density,  given  by 

n  b„ 

(?) 

°  j 

1  That  is,  if  AapY  are  the  components  of  a  rank-three  tensor,  then  per (A)apy  =  AapY  +  Aycp  +  Apya. 
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and  the  momentum  density,  given  by 

«  b„ 

=  m“c«J  N<‘.  i  ■  (8) 

«  i 

We  do  not  consider  a  hydrodynamic  equation  for  energy  since,  in  the  incompressible  limit,  that  decouples  from 
the  evolution  equations  for  mass  and  momentum.  That  is,  the  incompressibility  condition  and  momentum  equation 
may  be  solved  for  p  and  u  without  need  for  the  energy  equation  though,  of  course,  the  reverse  is  not  true.  In  other 
words,  the  energy  becomes  a  passive  variable  with  respect  to  the  dynamics  in  the  incompressible  limit,  so  we  need 
not  consider  it  here. 

In  keeping  with  earlier  work  [1],  we  assume  that  the  H  function  is  of  trace  form 
ba 

(9) 

a  i=l 

where  h'(x)  >  0  for  x  >  0.  If  we  extremize  this  under  the  assumption  that  p  and  pu  are  conserved,  we  find  the 
equilibrium  distribution  function 


Naj  =  <t>(P-ma  +  p  ■  macaj), 


where  p  and  p  are  the  Lagrange  multipliers  determined  by  the  constraints,  and  <p  the  inverse  function  of  h'. 

We  expand  the  equilibrium  distribution  to  second  order  in  Mach  number  formally,  using  p  as  an  expansion 
parameter 


NaSf  =  00 ima)  +  mac//(pma)P  ■  C aJ  +  \m2a<j)" (pma)PP  :  c ajcaj. 


Eqs.  (7)  and  (8)  are  then  used  to  derive  the  constraints 


p  =  ^ maNaj  =^'Ylma<p{pma)bQia  +  ^ PP  ■  ^ w^"(p,mn)b2,a, 

a,j  a  a 

pu  =  maCajN{aej  =  p  2, a- 


Under  the  assumption  that  bi,a  —  bixUli  for  l  <  4,  the  above  equations  become 


P  =  Yl,mab^a(j)(pma)  +  y  '22mlb2,a<P"(P'ma), 

a  a 

pu  =  PlLj  mlb2,a<t>' (t4ma). 


The  second  of  these  yields 


mlb2,a<P'(pma)  ’ 


and  the  first  then  yields 


V"1  ..  ,, .  s  ,  A2  'Eamlb2,a4>"(^ma) 

P —  /  1mab0%a4>(P'ma)  0  oj,  AJt  wl ' 

v  .  2  lEamib2,a4>(P-ma)V 


(17) 
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To  proceed,  we  define  the  functions 

®n(z)  =  ^mabn;a<Kzma)- 

Cl 

Note  that  the  f  th  derivative  of  these  is  given  by 
®n}  (z)  =  rnla+ 1  bn,a<t>(l)  C zma )■ 


In  terms  of  these,  p  may  be  written 

p2u?  0j(p)  ,  p2n2  0j(p) 

'  “  + + »lr^tW 


=  4>o  (k 


P2U2  Q^ip) 

2  <P{»(/*)[<^(#*)]3 


The  Lagrange  multiplier  p  is  then  obtained  by  first  writing 

p,  =  <Pq\p)  -  P~  TT r^l~y',2 ’  (21) 

and  applying  this  equation  to  itself  iteratively,  until  the  right-hand  side  no  longer  contains  p  explicitly.  We  find 

_ fx&'m. _ .  (22) 

It  follows  that  the  first  term  on  the  right  in  Eq.  (11)  may  be  written 

,2  v  ,  \  P2u20'2(‘Pol(p))ma<l>,(<Pol(p)ma ) 

<p(pma)  =  <£(<P0  (p)ma)  - - -_i~  ,^-1,  ,,l2  '  (23) 

Finally,  the  Lagrange  multiplier  ft  is  given  by 

(24) 

0'2(P)  .  <P'2(<Pq\p)) 

and  we  are  ready  to  write  down  the  complete  Mach-expanded  equilibrium  distribution  function,  expressed  in  terms 
of  the  conserved  quantities  p  and  pu 

Ar(eq)  ,2^,-1 2  s  ,  ,  ma<P'(0o\p)ma) 

NaJ  =  <H0O  (p)ma)  +  _s_T^r.u  .  caJ 

+  2[4($“;))P  :  .  (25) 


4.  Hydrodynamic  equations 


We  now  insert  the  form  of  the  equilibrium  distribution  derived  in  Section  3  into  the  BGK  kinetic  equation,  Eq.  (1), 
and  perform  the  Chapman-Enskog  asymptotic  analysis  for  small  Knudsen  number.  In  fact,  we  take  the  Knudsen 
and  Mach  numbers  to  be  of  the  same  order,  as  is  appropriate  for  the  incompressible  Navier-Stokes  equations.  The 
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details  (if  this  analysis  arc  presented  in  Appendix  A.  The  resulting  hydrodynamic  equations  are  V  •  u  =  0  and  the 
incompressible  Navier-Stokes  equation 


3u 

Yt 


+  gu  Vu  =  —V  P  +  vV2u, 


(26) 


where  we  have  defined  the  scalar  pressure 


the  kinematic  viscosity 


and  the  factor  multiplying  the  convective  derivative 
g  (^)2‘ 


P  =  <P2  + 


<P0<P'i  <P0®2 

<P'0' 


(27) 


(28) 


(29) 


Here,  all  of  the  functions  0n  are  understood  to  be  evaluated  at  1  (p).  We  note  that  the  correct  form  of  the  convective 

derivative,  and  therefore  Galilean  invariance,  is  recovered  when  g  =  1. 

It  should  be  noted  that  the  expression  for  i>  above  is  specific  to  the  choice  of  a  BGK  collision  operator  with 
constant  r,  but  the  results  for  g  and  P  are  more  general  than  that.  Once  we  determine  the  form  for  the  function  h 
that  will  make  g  —  1,  we  could  use  it  to  construct  a  variable-r  entropic  lattice  Boltzmann  collision  operator,  with 
relaxation  parameter  k  as  described  above,  and  perform  a  Chapman-Enskog  analysis  of  that.  The  important  result 
that  g  —  1  would  be  unaltered,  as  would  the  equation  of  state,  but  the  expression  for  the  viscosity  would  be  different, 
approaching  zero  as  k  approaches  unity.  We  shall  not  provide  the  details  of  this  entropic  collision  operator  in  this 
paper,  but  rather  focus  solely  on  the  determination  of  h. 


5.  Examples 


The  requirement  of  Galilean  invariance  is  then  g  =  1,  or 

0Q04  _ 

(0'2)2  ■ 


(30) 


In  this  section,  we  solve  this  equation  for  four  different  example  lattice  models.  The  first  is  the  single-speed  lattice 
Boltzmann  equation  with  a  Bravais  lattice;  here  we  make  contact  with  previous  results.  In  our  second  example,  we 
generalize  this  to  include  a  zero- velocity  component  to  the  distribution — the  case  of  so-called  “rest  particles” — with 
the  same  mass  as  the  moving  particles.  Here  we  show  that  the  solution  for  h  still  has  power-law  form,  and  present 
an  analytic  form  for  it  that  generalizes  the  result  for  a  Bravais  lattice.  In  our  third  example,  we  allow  the  rest  particle 
to  have  a  different  mass  from  the  moving  particles.  Remarkably,  we  can  show  that  a  power-law  solution  for  h  is  still 
obtained,  that  the  power  required  satisfies  a  certain  transcendental  equation,  and  that  this  equation  is  guaranteed  to 
have  a  solution  for  that  power.  Finally,  we  treat  the  general  case,  subject  to  the  restrictions  on  our  choice  of  lattice 
that  were  described  above.  We  show  that  the  power-law  form  again  holds,  and  that  the  power  required  again  satisfies 
a  certain  transcendental  equation;  in  this  general  case,  however,  we  are  unable  to  demonstrate  the  existence  of  a 
solution  for  that  power. 
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5.  /.  Single-speed  models 

In  the  case  of  a  single-speed  model,  the  index  a  takes  on  only  one  value,  which  we  may  suppress.  We  have 

0o(z)  =  mbo4>(.zm), 

bi 

02  (z)  =  mb2<t>(.zm)  —  ■rL0o(z), 

bo 

b4 

@4  (z)  =  mb4<j>(zm)  =  —  0q{z), 
bo 

so  Eq.  (30)  yields 

<fi(x)<p"(x)  =  X2[<fi'(x)]2, 

where  we  have  defined  the  new  independent  variable  x  =  zm,  and  the  coefficient 

X-  -2 
~Jbob4 

Eq.  (34)  is  a  generalization  of  the  differential  equation  derived  in  [1].  This  equation  has  solution 
<f>(x)  =  A(x~  B)1/(1~k2\ 

where  A  and  B  are  arbitrary  constants.  The  inverse  function  of  this  must  then  be 


/  7  \  1  —A2 

which  integrates  to  yield 


h(z)  = 


ho  +  Bz  +  A  In  z 


ho  +Bz 


if  X  =  72, 


otherwise. 


The  leading  linear  function  of  z  is  unimportant,  as  it  results  in  nothing  more  than  a  constant  addition  to  the  H 
function;  likewise,  the  multiplicative  constant  in  front  of  the  remaining  piece  is  unimportant  and  may  be  set  to  unity. 
What  remains  has  the  form  of  a  Tsallis  entropy  with  q  =  2  —  X2  for  A.  ^  and  a  Burg  entropy  otherwise.  For  a 
Bravais  lattice,  we  have  X  =  VI  +  2/D,  so  that  q—  1  —  2/D,  as  reported  in  earlier  work  [1]. 

5.2.  Models  with  rest  particles  of  the  same  mass 


To  understand  the  complications  that  arise  when  more  than  one  speed  is  used  in  entropic  lattice  Boltzmann 
models,  we  next  consider  the  addition  of  br  “rest  particles”  to  the  bm  single-speed  model  described  above.  Here 
we  have  introduced  the  speed  labels  a  —  r  for  the  rest  particles,  and  a  =  m  for  the  moving  particles.  We  note  that 
bo ,r  =  br,  but  that  &2,r  =  b4j  =  0  since  the  rest  particles  have  zero  speed.  It  follows  that 

<P0(z)  =  rnmbo,m<t>(zmm)  +  mrbo,r4>(zmr),  (39) 

(z)  =  rnmb2,mf(zmm),  (40) 

04  {z)  =  rnmb4<m(t>{zmm),  (41) 
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so  Eq.  (30)  yields 

[(p{x)  +  y4>(nx)]<p"(x)  =  X2[(j>'(x)f, 

where  we  have  defined  the  new  independent  variable  x  —  zm,  and  the  coefficients 
X  =  bxm 

■Jbo 


y,^.  (45, 

Wlm  O0,m 

If  the  rest  particles  have  the  same  mass  as  the  moving  particles,  then  p  —  1 ,  and  Eq.  (42)  reduces  to 
A2 

4>{x)<t>"{x)  =  -~[4>\x)f.  (46) 

1  +  Y 

This  is  of  the  same  form  as  the  equation  obtained  with  no  rest  particles,  except  with  the  substitution  X2  ->  X2/(l  +  y). 
It  follows  that  the  required  H  function  again  has  the  form  of  the  Tsallis  entropy  with 

X2 

4  =  2-——  (47) 

1  +  Y 

with  the  Burg  entropy  recovered  in  the  special  case  that  X2  =  2(1  +  y). 

For  a  Bravais  lattice,  the  expression  for  q  reduces  to 

,  =  (48) 

When  the  number  of  rest  particles,  &o,r  goes  to  zero,  this  reduces  to  the  single-speed  result. 

5.3.  Models  with  rest  particles  of  a  different  mass 


If  the  moving  particle  mass  and  rest  particle  mass  differ,  then  1,  and  Eq.  (42)  is  an  example  of  a  functional 
differential  equation  [13].  Such  equations  relate  the  dependent  variable  and  its  derivatives  at  more  than  one  value 
of  the  independent  variable.2  In  this  case,  the  value  of  the  dependent  variable  <j>  and  those  of  its  derivatives  at  x  are 
related  to  its  value  at  /ix.  Such  equations  are  difficult  to  solve,  even  when  they  are  linear  and  the  relevant  values 
of  the  independent  variable  are  related  by  additive  shifts;  our  equation  is  nonlinear,  and  the  relevant  values  of  the 
independent  variable  are  related  by  a  multiplicative  shift.  Remarkably,  in  spite  of  all  this,  our  equation  admits  an 
analytic  solution,  namely  the  power-law 

4>(x)  —  Ax^,  (49) 


where  A  is  an  arbitrary  constant,  and  f  solves  the  transcendental  equation 

1  +  YPfi  _  f 
X2  ~  p  -  1 ' 


(50) 


2  Here  we  are  using  the  term  functional  differential  equation  in  a  sense  most  often  used  by  mathematicians  [13],  Physicists  often  use  the  term 
to  refer  to  differential  equations  that  involve  what  they  call  functional  derivatives ,  which  are  what  mathematicians  call  Frechet  derivatives.  Lest 
there  be  any  confusion,  that  is  most  emphatically  not  what  we  are  talking  about  here. 
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Tlic  inverse  function  of  (p |S 

and  this  integrates  to  give 


Ap  /  z  \  1+1//3 

ft(z)  =  Ao  +  TJL(_)  , 


The  contribution  to  H  from  a  given  site  is  then 


n  b„  (  n  ba  \  i  /  a  \  «  ba 

E  E  = *o  ^E  E  >)  +  jro  (r^)  E  E  <r- 


Since  A  and  ho  are  arbitrary,  it  is  clear  that  this  can  be  brought  into  correspondence  with  the  Tsallis  entropic  form, 
which  is  a  linear  function  of  the  sum  over  N%  ..  We  thus  identify  q  =  1  +  1  //?,  where  p  must  be  found  by  solving 
the  transcendental  equation,  Eq.  (50). 

Since  we  would  like  the  distribution  function  <p  to  decrease  with  increasing  argument,  we  are  interested  in  solutions 
with  p  <  0.  If  X  >  1  and  /z  >  1,  such  a  solution  will  always  exist.  To  see  this,  note  that  the  function  k2p/(p  —  1) 
is  zero  when  p  =  0  and  approaches  A2  >  1  as  p  -*■  — oo;  then  note  that  the  function  1  +  is  1  +  y  >  1  when 
P  —  0  and  approaches  1  as  p  -»  — oo.  Since  both  curves  are  continuous,  they  must  cross  for  some  p  <  0. 

5.4.  The  general  case 

Encouraged  by  the  last  example,  we  may  check  to  see  if  a  power-law  solution  always  exists  for  <p.  If  we  assume 
4>(x)  =  Ax^,  then 

(z)  =  mabn,aMzma)P  =  Azp  (54) 

a  a 

Inserting  this  into  Eq.  (30),  we  quickly  find  that  must  satisfy  the  transcendental  equation 


(Eqml+%,a)(Ea'nl+%,a)  p 


The  existence  of  power-law  solutions  for  <p  then  hinges  on  the  existence  of  solutions  to  this  transcendental  equation 
with  p  <  0.  Though  we  saw  that  such  solutions  exist  for  single-speed  models  and  models  with  a  single-speed 
plus  rest  particles,  it  seems  difficult  to  draw  conclusions  about  the  existence  of  more  general  solutions  to  this 
transcendental  equation,  and  we  leave  the  matter  to  future  inquiry. 


6.  Conclusions 

Previous  work  [1]  has  established  that  the  requirement  of  Galilean  invariance  determines  the  form  of  the  H 
function  for  single-speed,  single-mass  entropic  lattice  Boltzmann  models  on  a  Bravais  lattice.  This  function  was 
found  to  have  the  form  of  a  Tsallis  entropy,  with  q  =  1—2 /D,  where  D  is  the  spatial  dimension.  In  this  study, 
we  have  extended  the  validity  of  this  result  to  models  with  multiple  speeds,  and  particles  of  different  masses.  We 
carried  out  the  full  Boltzmann/Chapman-Enskog  analysis  for  such  models,  and  applied  our  results  to  four  examples. 
The  first  was  the  single-speed  model,  to  verify  that  we  could  reproduce  the  earlier  result.  The  second  was  the  same 
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model  with  the  addition  of  rest  particles  of  the  same  mass;  here  we  showed  that  the  only  effect  was  to  change  the 
value  of  q.  Our  third  and  fourth  examples  involved  particles  of  more  than  one  mass;  here  we  showed  that  h  must 
satisfy  a  certain  functional  differential  equation,  that  this  equation  has  a  solution  in  power-law  form,  and  that  the 
power  is  determined  by  a  transcendental  equation.  Thus,  we  have  shown  that  the  power-law  form  for  h  is  extremely 
robust. 

We  note  that  our  choice  of  the  form  of  H  differs  from  that  of  Karlin  et  al.  [6],  which  is  of  the  form  H  = 
Yjj  Nj  \n(Nj/Wi),  where  the  W;  are  speed-dependent  weights.  These  weights  are  equal  to  the  global  equilibrium 
at  zero  flow;  thus,  when  Nj  =  Wj  they  have  H  =  0.  Thus,  by  allowing  weighted  contributions  to  H,  they  found 
solutions  for  which  h  has  the  form  of  a  (relative)  Boltzmann  entropy;  by  contrast,  the  present  work  assumes  uniform 
contributions  to  H,  and  finds  solutions  for  which  h  is  not  a  Boltzmann  entropy.  Both  approaches  are  capable  of 
yielding  Galilean-invariant  hydrodynamics.  A  more  general  form  for  H  that  will  subsume  both  approaches  as  special 
cases,  remains  an  interesting  theoretical  challenge. 

Another  outstanding  problem  involves  the  restriction  that  we  imposed  on  the  choice  of  lattice.  Though  the 
requirement  that  the  second-rank  and  fourth-rank  tensors  formed  from  sums  of  outer  products  of  the  velocities 
be  isotropic  within  each  speed  class  separately  allowed  for  a  simplified  analysis,  it  seems  unnecessary.  There  are 
known  lattice  Boltzmann  models  for  which  this  is  not  true,  but  which  nevertheless  yield  the  isotropic  Navier-Stokes 
equations  in  the  hydrodynamic  limit;  this  is  because  the  union  of  velocity  vectors  across  all  speed  classes  does 
satisfy  the  isotropy  requirements.  For  example,  the  very  popular  D2Q9,  D3Q15  and  D3Q19  models  fall  into  this 
category  [3], 

The  analysis  of  the  transcendental  equation  for  ft,  Eq.  (55),  to  see  under  what  circumstance  there  exist  solutions 
with  ft  <  0,  remains  yet  another  outstanding  problem. 

Finally,  in  addition  to  these  technical  challenges,  it  would  be  useful  and  enlightening  to  have  some  physical 
reason  for  the  appearance  of  the  Tsallis  entropic  form  in  this  context.  This  form  has  often  been  reported  as  arising 
from  a  lack  of  ergodicity,  or  a  fractal  spatiotemporal  structure.  There  is  no  clear  reason  to  believe  that  either  of 
these  ingredients  are  present  in  the  current  context;  yet  the  Tsallis  entropic  form  appeared  quite  naturally  from  our 
mathematical  development.  Thus,  a  clear  and  illuminating  physical  interpretation  of  our  result  remains  perhaps  the 
most  important  outstanding  challenge  that  we  leave  for  future  work. 
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Appendix  A.  Lattice  Boltzmann  equation  and  Chapman-Enskog  analysis 

Using  the  equilibrium  distribution  function  derived  in  Section  3,  the  kinetic  equation,  Eq.  (1),  may  be  rewritten 
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+  r  exp(c„.j  •  V)  exp  ^A/^  _  1  J  jV«.,/(x*  0  =  Naj(x>  ?)> 


which  has  the  formal  solution 


exp(cflj  •  V)  exp  ~  1  j  Naf(x’  ‘ 


Naj (x,  t)  -  1  +  r 


We  introduce  the  order  parameter  e  by  the  following  prescription,  appropriate  for  viscous  incompressible  flow 
(Mach  and  Knudsen  numbers  of  order  e) 


A  t  -»  ezA  t, 

and  we  allow  for  the  equilibrium  distribution  to  be  ordered  in  Mach  number 


AfCeq)  =  YV/V(cq’f) 

a,j  L-j  a,j  ' 


The  result 


Wa  J  =  1  +  r 


exp(ecfl>7-  -  V)  exp  (e2At^  -  1 J  J 


may  be  solved  perturbatively  by  ordering  Naj  in  the  expansion  parameter  e 
oo 

^=X>‘0 


At  orders  zero,  one  and  two,  we  find 
a.)  ~  a,j  ’ 

Jyaj  ~  lyaJ  TC“-J  ylya,j  ’ 


Kj  =  Naj2)  -  tcaj  •  -  r  [a*|  -  (r  -  I)  caJcaJ  :  w]  Ag'0). 


(A.  10) 


(All) 


We  shall  use  these  to  derive  the  leading  gradient  corrections  to  the  distribution  function,  and  insert  these  into  the 
conservation  equations  to  arrive  at  the  hydrodynamical  equations  for  the  system. 

We  now  proceed  to  the  Chapman-Enskog  analysis.  The  ordering  of  our  local  equilibrium  distribution  function  is 


Kj  =  W  l(P)ma), 

„(eq,l)  rna4>’(<pQl(p)ma) 

~~  1/  vs  Pu'caJ< 


^(eq,2)  _ 


<P'2^ol(p)) 

2  w™))]2 : 


(A. 12) 


(A.13) 


(A. 14) 


(A.  15) 


(A.  16) 


(A.  17) 
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Inserting  Eqs.  (A.  1 2)-(A.  14)  into  Eqs.  ( A.9)-(A.  I  1 ),  we  get 

^0)j=<P(<Po'(p)mu),  (A.  15) 

,,(l)  ma<p'(<pQ\p)ma) 

w <A  6) 

N“')j  =  2[4>^(4>o1(p))]2  1  ^  (0°  1(p)ma)CajCa'j  ~  ma4>  (<P°l{p)ma)  cp^-'ip))1 

nna(t>'(<Pol(p)ma ) 

- ^TT^T -—-PCaJCaJ  :  (Vu).  (A.  17) 

0>)) 

The  solution  to  the  kinetic  equation  to  second  order  is  then  the  sum  of  these 

Na,j  =  Nj®  +  eN™  +  €2N^j  +  0(€3),  (A.  18) 

and  the  conservation  equations  are  obtained  by  taking  moments  of  the  Taylor-expanded  kinetic  equation  to  second 
order 

+  C aJ  ■  V Naj  +  IrCajCaj  :  VVNaJ  =  £2aJ,  (A.  19) 

01  z 

where  S2a,j  denotes  the  collision  operator.  The  mass  moment  yields 

~  +  V  -  (pn)  +  \ VV  :  ^  maCa.jCa,jNa,j^  =  0,  (A.20) 

while  the  momentum  moment  yields 


(A.  18) 


(A.19) 


(A.20) 


ainc  me  liiuineiiLuiii  niuinciu  yieiub 

d  (  n  bo  \  i  /  n  b„  \ 

—  (pu)  +  V  •  I  ^2Y^macajCajNajJ  +  -VV  :  I  EE  maCa,jCa,jca,jNaj  I  —  0.  (A.21) 

Ve  now  evaluate  the  second  and  third  moments  which  appear  in  these  equations,  order  by  order.  For  the  second 
- - - £2 _ 1 


noments,  we  find 

f^J^rriaCajCajN^j  =  02, 


T.  y  maca,jcajNaj  —  0, 


(2)  [&0&A  0O&U  pu2  K  2  T 

£X>«c ajcaJNaJ=  — 1+  —  puu-r^p[(V.u)l+Vu+(Vu)  ], 

Cl  J  “2  v  a  "*  2 


(A.  22) 


(A.23) 


(A.24) 


lere  the  superscript  T  denotes  “transpose”,  and  the  arguments  of  all 
ird  moments,  we  find 


the  0’s  are  understood  to  be  0O  1  (p).  For  the 


n  ba 

T.  y,  maCa,jCa,jCa,jNaJ  =  ®> 


(A.25) 
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0' 

J2J2m‘<C“JC“-jC“JNaJ  =  ®  u)> 

a  j  2 


(A. 26) 


n  b„ 

J2  £  ”'ac(,jCa.jCa.jN™  =  0,  (A. 27) 

“  J 

where,  as  noted  earlier,  “per”  denotes  a  sum  over  all  symmetric  permutations  of  tensor  components. 

We  now  insert  the  above  results  into  Eqs.  (A.20)  and  (A.21),  and  we  adhere  to  the  incompressible  limit  so  that 
time  and  space  derivatives  of  p  are  of  higher  order  in  e  and  hence  ignored.  The  hydrodynamic  equations  given  in 
Section  4  follow  immediately. 
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